CC=================================================================CC
CC                                                                 CC
CC  Subroutine IRSRFFT(X,M):                                       CC
CC      A inverse real-valued, in-place, split-radix FFT program   CC
CC      Decimation-in-frequency, cos/sin in second loop            CC
CC      and computed recursively                                   CC
CC      Symmetric input in order:                                  CC
CC              [ Re(0),Re(1),....,Re(N/2),Im(N/2-1),...Im(1)]     CC
CC      The output is real-valued                                  CC
CC                                                                 CC
CC  Input/output                                                   CC
CC      X    Array of input/output (length >= N)                   CC
CC      M    Transform length is N=2**M                            CC
CC                                                                 CC
CC  Calls:                                                         CC
CC      IRSTAGE, RBITREV                                           CC
CC                                                                 CC
CC  Author:                                                        CC
CC      H.V. Sorensen,   University of Pennsylvania,  Oct. 1985    CC
CC                       Arpa address: hvs@ee.upenn.edu            CC
CC  Modified:                                                      CC
CC      F. Bonzanigo,    ETH-Zurich,                  Sep. 1986    CC
CC      H.V. Sorensen,   University of Pennsylvania,  Mar. 1987    CC
CC                                                                 CC
CC  Reference:                                                     CC
CC      Sorensen, Jones, Heideman, Burrus :"Real-valued fast       CC
CC      Fourier transform algorithms", IEEE Tran. ASSP,            CC
CC      Vol. ASSP-35, No. 6, pp. 849-864, June 1987                CC
CC      Mitra&Kaiser: "Digital Signal Processing Handbook, Chap.   CC
CC      8, page 491-610, John Wiley&Sons, 1993                     CC
CC                                                                 CC
CC      This program may be used and distributed freely provided   CC
CC      this header is included and left intact                    CC
CC                                                                 CC
CC=================================================================CC
        SUBROUTINE IRSRFFT(X,M)
        REAL X(1)
C-------L shaped butterflies----------------------------------------C
        N = 2**M
        N2 = 2*N
        DO  10 K = 1, M-1
            N2 = N2/2
            N4 = N2/4
            CALL IRSTAGE(N,N2,N4,X(1),X(N4+1),X(2*N4+1),X(3*N4+1))
 10     CONTINUE
C-------Length two butterflies--------------------------------------C
        IS = 1
        ID = 4
 70         DO  60 I1 = IS,N,ID
                T1      = X(I1)
                X(I1)   = T1 + X(I1+1)
                X(I1+1) = T1 - X(I1+1)
 60         CONTINUE
            IS = 2*ID - 1
            ID = 4*ID
        IF (IS.LT.N) GOTO 70
C-------Digit reverse counter---------------------------------------C
        CALL RBITREV(X,M)
C-------Divide by N-------------------------------------------------C
        DO  99 I=1,N
            X(I) = X(I)/N
 99     CONTINUE
        RETURN
        END
CC=================================================================CC
CC                                                                 CC  
CC  Subroutine IRSTAGE - the work-horse of the IRFFT               CC
CC       Computes a stage of an inverse real-valued split-radix    CC
CC       length N transform.                                       CC
CC  Author                                                         CC
CC       H.V. Sorensen,   University of Pennsylvania,  Mar. 1987   CC
CC                                                                 CC
CC      This program may be used and distributed freely provided   CC
CC      this header is included and left intact                    CC
CC                                                                 CC
CC=================================================================CC
        SUBROUTINE  IRSTAGE(N,N2,N4,X1,X2,X3,X4)
        DIMENSION X1(1),X2(1),X3(1),X4(1)
        N8 = N4/2
        IS = 0
        ID = 2*N2
 10         DO  20  I1 = IS+1,N,ID
                T1     = X1(I1) - X3(I1)
                X1(I1) = X1(I1) + X3(I1)
                X2(I1) = 2*X2(I1)
                T2     = 2*X4(I1)
                X4(I1) = T1 + T2
                X3(I1) = T1 - T2
 20         CONTINUE
            IS = 2*ID - N2
            ID = 4*ID
        IF (IS .LT. N) GOTO 10
C
        IF  (N4-1) 100,100,30
 30     IS = 0
        ID = 2*N2
 40         DO  50 I1 = IS+1+N8,N,ID
                T1    = (X2(I1) - X1(I1))*1.4142135623730950488
                T2    = (X4(I1) + X3(I1))*1.4142135623730950488
                X1(I1) = X1(I1) + X2(I1)
                X2(I1) = X4(I1) - X3(I1)
                X3(I1) = -T2-T1
                X4(I1) = -T2+T1
 50         CONTINUE
            IS = 2*ID - N2
            ID = 4*ID
        IF (IS .LT. N-1) GOTO 40
C
        IF  (N8-1) 100,100,60
 60     E  = 6.283185307179586/N2
        SS1 = SIN(E)
        SD1 = SS1
        SD3 = 3.*SD1-4.*SD1**3
        SS3 = SD3
        CC1 = COS(E)
        CD1 = CC1
        CD3 = 4.*CD1**3-3.*CD1
        CC3 = CD3
        DO  90  J = 2,N8
            IS = 0
            ID = 2*N2
            JN = N4 - 2*J + 2
 70             DO  80 I1=IS+J,N,ID
                    I2 = I1 + JN
                    T1     = X1(I1) - X2(I2)
                    X1(I1) = X1(I1) + X2(I2)
                    T2     = X1(I2) - X2(I1)
                    X1(I2) = X2(I1) + X1(I2)
                    T3     = X4(I2) + X3(I1)
                    X2(I2) = X4(I2) - X3(I1)
                    T4     = X4(I1) + X3(I2)
                    X2(I1) = X4(I1) - X3(I2)
                    T5 = T1 - T4
                    T1 = T1 + T4
                    T4 = T2 - T3
                    T2 = T2 + T3
                    X3(I1) =  T5*CC1 + T4*SS1
                    X3(I2) = -T4*CC1 + T5*SS1
                    X4(I1) =  T1*CC3 - T2*SS3
                    X4(I2) =  T2*CC3 + T1*SS3
 80             CONTINUE
                IS = 2*ID - N2
                ID = 4*ID
            IF (IS .LT. N) GOTO 70
C
            T1  = CC1*CD1 - SS1*SD1
            SS1 = CC1*SD1 + SS1*CD1
            CC1 = T1
            T3  = CC3*CD3 - SS3*SD3
            SS3 = CC3*SD3 + SS3*CD3
            CC3 = T3
 90     CONTINUE
 100    RETURN
        END
CC=================================================================CC
CC                                                                 CC
CC Subroutine RBITREV(X,M):                                        CC
CC      Bitreverses the array X of length 2**M. It generates a     CC
CC      table ITAB (minimum length is SQRT(2**M) if M is even      CC
CC      or SQRT(2*2**M) if M is odd). ITAB need only be generated  CC
CC      once for a given transform length.                         CC
CC                                                                 CC
CC Author:                                                         CC
CC      H.V. Sorensen,   University of Pennsylvania,  Aug. 1987    CC
CC                       Arpa address: hvs@ee.upenn.edu            CC
CC References:                                                     CC
CC      D. Evans, Tran. ASSP, Vol. ASSP-35, No. 8, pp 1120-1125    CC
CC      Aug. 1987                                                  CC
CC                                                                 CC
CC      This program may be used and distributed freely provided   CC
CC      this header is included and left intact                    CC
CC                                                                 CC
CC=================================================================CC
        SUBROUTINE RBITREV(X,M)
        DIMENSION X(1),ITAB(256)
C-------Initialization of ITAB array--------------------------------C
        M2 = M/2
        NBIT = 2**M2
        IF (2*M2.NE.M) M2 = M2 + 1
        ITAB(1) = 0
        ITAB(2) = 1
        IMAX = 1
        DO  10 LBSS = 2, M2
            IMAX = 2 * IMAX
            DO  10 I = 1, IMAX
                ITAB(I)      = 2 * ITAB(I)
                ITAB(I+IMAX) = 1 + ITAB(I)
 10     CONTINUE
C-------The actual bitreversal--------------------------------------C
        DO  20 K = 2, NBIT
            J0 = NBIT * ITAB(K) + 1
            I = K
            J = J0
            DO  20 L = 2, ITAB(K)+1
                T1   = X(I)
                X(I) = X(J)
                X(J) = T1
                I = I + NBIT
                J = J0 + ITAB(L)
 20     CONTINUE
        RETURN
        END
